function [ results_fnames ] = osl_run_first_level_continuous( oat )

% [ results_fnames ] = osl_run_first_level_continuous( oat )
%
% takes in an OAT, which needs to be setup by calling oat=osl_setup_oat(S), struct
% and runs first level continuous time GLM or seed based connectivity
%
% This function should normally be called using osl_run_oat(oat);
%
% MWW 2012

global OSLDIR;

use_class_probs=0;

first_level=oat.first_level;
source_recon_name='source_recon';

% check if user has asked not to do GLM
if(~isfield(first_level,'doGLM')),
    first_level.doGLM = 1;
end

if(first_level.time_average)
    error('Time averaging over whole window is not appropriate on continuous data');   
end;

if ~isempty(first_level.tf_downsample_factor)
    first_level.tf_downsample_factor=1;
    warning('Turning off downsampling - not implemented');
end;

% set the downsample factor.  If doing a tf, read tf_downsample_factor.  If
% doing time domain, read time_downsample_factor
switch first_level.tf_method
    case 'none' % use time_downsample_factor
        error('Continuous data time-wise GLM only handles time-freq analysis');
        
        if ~isempty(first_level.time_downsample_factor)
            ds_factor = 1./first_level.time_downsample_factor;
        else
            ds_factor = [];
        end
    otherwise  % use tf_downsample_factor
        if ~isempty(first_level.tf_downsample_factor)
            ds_factor = 1./first_level.tf_downsample_factor;
        else
            ds_factor = [];
        end
end % switch first_level.tf_method

if(~isfield(first_level,'connectivity_seed_mni_coord')),
    first_level.connectivity_seed_mni_coord=[];
end;

Ntrials   = 1;    
   
clear stats;

for subi_todo=1:length(first_level.sessions_to_do),   
        
    sub=first_level.sessions_to_do(subi_todo);
    
    %sub/length(first_level.submatfiles),
    
    %% load in the beamformer result for this session
    source_recon_results_fname = oat.source_recon.results_fnames{sub};
    source_recon_results=osl_load_oat_results(oat,source_recon_results_fname);
    
    % get D and update path in case OAT dir has been moved
    D=source_recon_results.BF.data.D;
    if(~strcmp([oat.source_recon.dirname],D.path)),        
        D=path(D,[oat.source_recon.dirname]);
        D.save;
        D=spm_eeg_load([D.path '/' D.fname]);
        osl_save_oat_results(oat,source_recon_results);    
    end;
    
    %% setup things that are common to all sessions
    if(subi_todo==1),
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up mask indices from mask/coords
        if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
            
            disp('Working in sensor space');
  
            %%%%%%%%%%%%%%%
            %% setup channels
            modality='MEG';
            chanindmeg = strmatch(modality, D.chantype); % need to even include bad channels in output so that it is the same for all sessions
       
            first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
                   

        else,
    
            if(isfield(first_level,'mni_coords'))    
                                
                current_level_mni_coord = first_level.mni_coords;
                
                % load in lower level coords
                lower_level_mni_coord=source_recon_results.mni_coord;               
                %lower_level_mask_fname  = [oat.source_recon.dirname '/' source_recon_name '_mask'];
                %[ lower_level_mni_coord, ~ ] = osl_mnimask2mnicoords(lower_level_mask_fname);

                [~, iA, mask_indices]=intersect(current_level_mni_coord, lower_level_mni_coord,'rows');

                % sort indices to be in the order of the lower level mni_coords:
                [~, gg]=sort(iA);
                indices_in_lower_level=mask_indices(gg);

                first_level_results.mask_indices_in_source_recon = indices_in_lower_level;
                first_level_results.mni_coord = lower_level_mni_coord(indices_in_lower_level,:);
                
                % the mni coordinates from the specified list
                str=['Desired MNI coordinates ' ];
                for vox=1:size(current_level_mni_coord,1),
                    str=[str ', [' num2str(current_level_mni_coord(vox,:)) ']'];
                end;
                disp(str);
                
                % the mni coordinates from the lower level mask
                str=['Nearest dipoles from lower level at MNI coordinates ' ];
                for vox=1:length(indices_in_lower_level),
                    indices_in_lower_level(vox)
                    str=[str ', [' num2str(lower_level_mni_coord(indices_in_lower_level(vox),:)) ']'];
                end;
                disp(str);
            
            else
                % setup std space brain
                first_level_results.gridstep=source_recon_results.gridstep;

                S=[];
                S.lower_level_mask_fname=[oat.source_recon.dirname '/' source_recon_name '_mask'];
                S.current_level_mask_fname=[oat.source_recon.dirname '/' oat.first_level.name '_mask'];
                S.current_level=first_level;
                S.lower_level_mni_coord=source_recon_results.mni_coord;
                S.lower_level_gridstep=source_recon_results.gridstep;

                [first_level_results.mask_indices_in_source_recon, first_level_results.mni_coord]=setup_mask_indices(S);
                lower_level_mni_coord=source_recon_results.mni_coord
                clear S;
            end;
        end;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup HMM 
        if strcmp(source_recon_results.recon_method,'none'),
            NK=1;
        else
            NK=numel(source_recon_results.BF.inverse.W.MEG); 
        end;
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up potential time frequency bases in tf_settings, note that any time or
        % frequency averaging will be done AFTER the time-frequency decomposition
        
        % set the frequency range
        freq_range = first_level.tf_freq_range; % range used for stat testing
        if isempty(freq_range)
            % if not set as a first_level option, default to the range used
            % for source_recon
            freq_range=source_recon_results.source_recon.freq_range;
            
            if isempty(freq_range)
                error('No frequency range has been specified. Set oat.source_recon.freq_range or oat.first_level.tf_freq_range');
            end % isempty(freq_range)
        end % if isempty(freq_range)
        
        source_recon_time_indices=find(source_recon_results.samples2use);
        source_recon_times = D.time(source_recon_time_indices);
        
        if isempty(first_level.time_range)
            first_level.time_range=source_recon_results.woi;
        end
        
        tf_settings = [];
        tf_settings.tf_method               = first_level.tf_method;
        tf_settings.tf_logtransform         = first_level.tf_logtransform;
        tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
        tf_settings.tf_freq_range           = freq_range; % set above
        tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
        tf_settings.raw_times               = source_recon_times;
        tf_settings.ds_factor               = 1; % do not downsample inside transform - instead we will do this later after the moving average low pass filtering
        tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
        tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
        tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
        tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;
        
        out = osl_tf_transform( tf_settings ); % passed without a data matrix, osl_tf_transform just gets container sizes
        
        if isfield(out, 'tf_morlet_basis')
            first_level.tf_morlet_basis = out.tf_morlet_basis;
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup time and frequency windows/indices for using later
          
        % setup result time dimension        
        first_level_time_indices  = ( out.tf_times >= first_level.time_range(1) ) & ( out.tf_times <= first_level.time_range(2) );
        first_level_results.times=1;       
        if first_level.time_average
            disp('Note time averaging irrelevant and so is ignored for continuous data time-wise GLM');
        end
        
        % setup result freq dimension
        first_level_results.frequencies = out.tf_freqs;
        if(first_level.freq_average),
            first_level_results.frequencies=mean(first_level_results.frequencies);
        end;
       
        % setup result spatial dimension
        Nvoxels_out=length(first_level_results.mask_indices_in_source_recon);

        clear out

        %%%%%%%%%%%%%%%%
        %% setup seed for seed-based connectivity analysis       
        if(~isempty(first_level.connectivity_seed_mni_coord)),
            
            if(first_level.tf_num_freqs>1)
                error('Continuous data seed-based connectivity only handles single freq bands at a time');
            end;

            mni_coords=lower_level_mni_coord;
            
            disp(['Using for seed: MNI coordinate ' num2str(first_level.connectivity_seed_mni_coord)]);
            
            % find the nearest index of the beamformed voxels to the specified mni_coord
            dists=(sqrt(sum((mni_coords-repmat(first_level.connectivity_seed_mni_coord,size(mni_coords,1),1)).^2,2)));
            
            [dist,first_level_results.seed_index]=min(dists);
            
            disp(['Nearest dipole for seed at MNI coordinate ' num2str(mni_coords(first_level_results.seed_index,:))]);                                              
                       
            first_level.contrast=[];
            first_level.contrast{1}=[1]';

        end;        
        
        %%%%%%%%%%%%%%%%
        %% results containers
        ntpts=1;
        nfreqs=length(first_level_results.frequencies);
        
        if(first_level.doGLM)
            first_level_results.stdcope=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs]);% num_voxels x  num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
            first_level_results.cope=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs]);% num_voxels x num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
            
            first_level_results.stdcope_by_state=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs,NK]);% num_voxels x  num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
            first_level_results.cope_by_state=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs,NK]);% num_voxels x num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
        end;

        first_level_results.pseudo_zstat_var=zeros([Nvoxels_out,1]);
        
    end; % end of if(sub==1)
    
    disp(['Reconstruct time courses and computing stats for dataset ' source_recon_results.fname]);    

    %%%%%%%%%%%%%%%
    %% setup channels
    modality='MEG';
    chanindmeg = strmatch(modality, D.chantype);
    chanind = setdiff(chanindmeg, D.badchannels);
    if isempty(chanind)
        error(['No good ' modality ' channels were found.']);
    end
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
        % for sensor space - also include bad sensors sensors 
        first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
        chanind = chanindmeg;
    end;    
    
    %%%%%%%%%%%%%%%
    %% load in sensor data
    first_level_results.D_sensor_data=D;
    first_level_results.chanind=chanind;
    
    triallist=1;
    
    sensor_data=D(chanind, source_recon_time_indices, 1);   
    class_prs=zeros(NK, length(source_recon_time_indices), length(triallist));

    classchanind=find(strcmp(D.chanlabels,'Class'));
    if(isempty(classchanind)),
        error(['No ''Class'' chanlabel in: ' D.fname]);
    else
        for si=1:NK,        
            class_samples_inds{si} = (D(classchanind,source_recon_time_indices, 1)==si);  

            state_tinds=(class_samples_inds{si}==1);                                                

            total_time=(source_recon_times(end)-source_recon_times(1))*size(state_tinds,3);
            state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
            txt=['State ' num2str(si) ' is active for ' num2str(state_time) 'secs'];

            if(state_time<10),
                warning([txt '. It will be excluded from the GLM analysis.']);
            else
                disp(txt);
            end;

            if(use_class_probs)
                classpr_chanind=find(strcmp(D.chanlabels,['ClassPr' num2str(si)]));
                class_prs(si,:,:) = D(classpr_chanind,source_recon_time_indices, triallist);
            end;            
        end; 
    end;

    
    %%%%%%%%%%%%%%%%
    %% setup the GLM design matrix
    % design matrix should be num_PEs x num_timepoints    
    if(~isempty(first_level.connectivity_seed_mni_coord)),

        % set the design matrix regressor to be the seed voxel power time course
        
        % recon seed time course
        S2=[];
        S2.source_recon_results=source_recon_results;
        S2.chanind=chanind;
        S2.class_samples_inds=class_samples_inds;
        S2.voxind=first_level_results.seed_index;
        S2.class_prs=class_prs;
        S2.use_class_probs=use_class_probs;
     
        S2.sensor_data=sensor_data;

        [ dat_seed wnorms_seed wnorm_seed wnorms_nai_seed wnorm_nai_seed] = get_voxel_recon_timecourse( S2 );
       
        dat_seed=dat_seed';
        wnorms_nai_seed=wnorms_nai_seed';
        
        dat_seed=dat_seed./sqrt(mean(squash(wnorms_nai_seed)));
        
        %% Transform the data - either downsampling in time domain, or a TF transform
            
        tf_settings = [];
        tf_settings.tf_method               = first_level.tf_method;
        tf_settings.tf_logtransform         = first_level.tf_logtransform;
        tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
        tf_settings.tf_freq_range           = freq_range; % set above
        tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
        tf_settings.raw_times               = source_recon_times;
        tf_settings.ds_factor               = 1; % do not downsample inside transform - instead we will do this later after the moving average low pass filtering
        tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
        tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
        tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
        tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;

        if isfield(first_level,'tf_morlet_basis')
            tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
        end

        out = osl_tf_transform( tf_settings , dat_seed' ); % do the transformation
       
        dat_seed_tf = permute(out.dattf,[2 3 1]);
                    
        first_level.design_matrix=dat_seed_tf';

        % add a regressor to model the mean and variance normalise
        % the seed power time course
        %processed_design_matrix=[(first_level.design_matrix); ones(size(first_level.design_matrix,2),1)'];  
        processed_design_matrix=[(first_level.design_matrix)];
              
        %first_level.contrast=[];
        %first_level.contrast{1}=[1 0]';
        
    else,
        
        if(~isfield(first_level,'design_matrix')),
                
            if(first_level.doGLM)
                error('Seed voxel or design matrix needs to be specified');
            else
                % setup dummy design matrix
                first_level.design_matrix=ones(1,length(source_recon_time_indices));
            end;
        end;
             
        processed_design_matrix=first_level.design_matrix(:,source_recon_time_indices);
    
    end;     
       
    if(size(processed_design_matrix)~=length(first_level_time_indices)),
        error(['Design matrix dimensions [npes x ntps] = ' num2str(size(first_level.design_matrix)) ' do not match the data: ntpts=' num2str(length(first_level_time_indices))]); 
    end;
    
    disp(['Design matrix dimensions [npes x ntps] = ' num2str(size(first_level.design_matrix))]);
        
    for pp=1:size(processed_design_matrix,1),
        processed_design_matrix2(pp,:)=processed_design_matrix(pp,first_level_time_indices);
    end;   
    processed_design_matrix=processed_design_matrix2;
                
    % apply moving average to design matrix
    if ~isempty(first_level.time_moving_av_win_size)
        fsample=1/(source_recon_times(2)-source_recon_times(1));        
        windowSize = first_level.time_moving_av_win_size*fsample;
                                
        for pp=1:size(processed_design_matrix,1),
            processed_design_matrix(pp,:)=moving(processed_design_matrix(pp,:),round(windowSize));
        end;        
    end;

    % temporally downsample                                     
    if ~isempty(first_level.tf_downsample_factor) && first_level.tf_downsample_factor~=1; 
        processed_design_matrix = osl_spm_resample(processed_design_matrix,ds_factor);        
    end;
    
    contrast_list=first_level.contrast;
    for c=1:length(contrast_list),
        contrast_list{c}=contrast_list{c}(:);
    end;

    %%%%%%%%%%%%%%%%
    %% do everything from now on one voxel at a time (for the sake of
    %% RAM) in the first level mask space                
    if(first_level.doGLM)
        disp('First level COPEs outputted will have dimension Nvoxels x Ntpts x Ncontrasts x Nfreqs:');
        tmp=[size(first_level_results.cope) 1 1];
        disp([num2str(tmp(1:4))]);
    end;
         
    %%%%%%%%%%%%%%%%
    %% setup design matrices that can be reused over voxels
    
    const_reg=0;        
    for p=1:size(processed_design_matrix,1),
        if(range(processed_design_matrix(p,:))==0),
            if(const_reg),
                warning('Multiple constant regressors in the design matrix!!!!?');
            end;
            const_reg=1;

        end;                                    
    end;

    if(const_reg && oat.first_level.do_glm_demean)
        warning('CAREFUL: are you sure you want the first_level.do_glm_demean flag on with constant regressors in the design matrix!!!!?');
    elseif(~const_reg && ~oat.first_level.do_glm_demean)
        warning('CAREFUL: are you sure you want no first_level.do_glm_demean flag on with no constant regressors in the design matrix!!!!?');
    end;       
        
    ntpts=1;
    x_stored=cell(NK,ntpts);
    pinvxtx_stored=cell(NK,ntpts);
    pinvx_stored=cell(NK,ntpts);
    
    for si=1:NK,
        state_tinds=(class_samples_inds{si}(first_level_time_indices)==1);

        total_time=(source_recon_times(end)-source_recon_times(1))*size(state_tinds,3);
        state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
 
        if(state_time>10),
            x=processed_design_matrix';
            x=x(permute(state_tinds(1,:),[2 1]),:);                        

            if(oat.first_level.do_glm_demean)
                for pp=1:size(x,2),
                    x(:,pp)=x(:,pp)-mean(x(:,pp));
                end;
            end;

            x_stored{si}=x;
            pinvxtx_stored{si}=pinv(x'*x);
            pinvx_stored{si}=pinv(x);
        end;
    end;
    
    contrast_list_mat=cell2mat(contrast_list);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% now loop over voxels
    ft_progress('init', 'etf');
    
    % debug:
    %[indfind, mni_coord]=nearest_vec(source_recon_results.mni_coord,[-32 -26 50]);
                            
    for indind=1:Nvoxels_out, % indexes brain space      
        
        ft_progress(indind/Nvoxels_out);                                   

        %voxind=first_level_results.mask_indices_in_source_recon(indind);
        
        S2=[];
        S2.source_recon_results=source_recon_results;
        S2.chanind=chanind;
        S2.class_samples_inds=class_samples_inds;
        S2.voxind=first_level_results.mask_indices_in_source_recon(indind);
        S2.sensor_data=sensor_data;
        S2.class_prs=class_prs;
        S2.use_class_probs=use_class_probs;

        [ dat wnorms wnorm wnorms_nai wnorm_nai] = get_voxel_recon_timecourse( S2 );

        %dat=dat'./sqrt(wnorms_nai');
        %wnorms=wnorms'./wnorms_nai';

        dat=dat';
        wnorms_nai=wnorms_nai';            
        dat=dat./sqrt(mean(squash(wnorms_nai)));

        first_level_results.pseudo_zstat_var(indind)=var(dat);

        if sum(isnan(squash(dat)))~=0,

            if(first_level.doGLM)
                first_level_results.cope(indind,:,:,:)=nan;
                first_level_results.stdcope(indind,:,:,:)=nan;         
                first_level_results.cope_by_state(indind,:,:,:,:)=nan;
                first_level_results.stdcope_by_state(indind,:,:,:,:)=nan; 
            end;
            warning('NANs in recon data');
        else
           
            if(~isempty(first_level.connectivity_seed_mni_coord)),
                if(first_level.connectivity_seed_regress_zerolag),

                    % orthoganalise the target voxel WRT the seed voxel (doing
                    % it this way round means that the TF transform and
                    % windowing does not need to be redone for the seed voxel)
                    if(NK>1),
                        
                        for i=1:NK,
                            tinds=(class_samples_inds{si}==1);
                            dat1=dat(tinds);
                            dat2=dat_seed(tinds);

                            dat(tinds)=dat1-dat2*(dat2\dat1);
                        end;     
                    else,
                        dat1=dat;
                        dat2=dat_seed;

                            
                        dat=dat1-dat2*(dat2\dat1);
                        
                    end;        
                    
                end;
            end;
            
            %% Transform the data - either downsampling in time domain, or a TF transform

            tf_settings = [];
            tf_settings.tf_method               = first_level.tf_method;
            tf_settings.tf_logtransform         = first_level.tf_logtransform;
            tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
            tf_settings.tf_freq_range           = freq_range; % set above
            tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
            tf_settings.raw_times               = source_recon_times;
            tf_settings.ds_factor               = 1; % do not downsample inside transform - instead we will do this later after the moving average low pass filtering
            tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
            tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
            tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
            tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;

            if isfield(first_level,'tf_morlet_basis')
                tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
            end

            out = osl_tf_transform( tf_settings , dat' ); % do the transformation

            dat_tf = permute(out.dattf,[2 3 1]);

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
            %% now do windowing over time 
            dat=dat_tf(first_level_time_indices,:);          
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
            %% now do averaging over time and frequency
            if(first_level.freq_average)
                dat=mean(dat,2);                
            end;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
            %% moving average within
            %% the window specified
            if ~isempty(first_level.time_moving_av_win_size)

                fsample=1/(source_recon_times(2)-source_recon_times(1));        
                windowSize = first_level.time_moving_av_win_size*fsample;

                dat=moving(dat,round(windowSize));             
            end

            % temporally downsample                                     
            if ~isempty(first_level.tf_downsample_factor) && first_level.tf_downsample_factor~=1;
                dat= osl_spm_resample(dat',ds_factor)';                  
            end;
            glm_input_times=linspace(first_level.time_range(1),first_level.time_range(2),length(dat));
                                                                        
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
            %% now do the GLM stats
            if first_level.doGLM
                                
                if(size(dat,1)>1)                       
                    
                    for si=1:NK,

                        state_tinds=(class_samples_inds{si}(first_level_time_indices)==1);
      
                        total_time=(source_recon_times(end)-source_recon_times(1))*size(state_tinds,3);
                        state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
                                              
                        if(state_time<10),
                            
                            first_level_results.cope_by_state(indind,:,:,:,si)=0;
                            first_level_results.stdcope_by_state(indind,:,:,:,si)=1e32;
                                    
                        else
                            
                            t=1;
                            x=x_stored{si};
                            pinvxtx=pinvxtx_stored{si};
                            pinvx=pinvx_stored{si};                               

                            for f = 1 : size(dat,2),

                                dattf=dat(state_tinds,t,f);

                                if(oat.first_level.do_glm_demean)
                                    dattf=demean(dattf);                                                        
                                end;

                                if(first_level.use_robust_glm),
                                    [b,stats2]=robustfit(x,dattf,'bisquare',4.685,'off'); % fits GLM robustly, Y=XB+E
                                    for c=1:length(contrast_list),

                                        if(strcmp(Sin.cope_type,'coape')),
                                            first_level_results.cope_by_state(indind,t,c,f,si)=contrast_list{c}'*abs(b);
                                        elseif(strcmp(Sin.cope_type,'cope')),
                                            first_level_results.cope_by_state(indind,t,c,f,si)=contrast_list{c}'*b;
                                        elseif(strcmp(Sin.cope_type,'acope')),
                                            first_level_results.cope_by_state(indind,t,c,f,si)=abs(contrast_list{c}'*b);
                                        else
                                            error('Invalid cope_type');
                                        end;
                                        first_level_results.stdcope_by_state(indind,t,c,f,si)=sqrt(contrast_list{c}'*stats2.covb*contrast_list{c});
                                    end;
                                    pe(t,:)=b; % num_timepoints x  num_pes
                                else
                                    
                                    
                                    [copeout, varcopeout, coapeout, dof,peout]=glm_fast_for_meg(dattf,x,pinvxtx,pinvx,contrast_list_mat,0);

                                    if(strcmp(oat.first_level.cope_type,'coape')),
                                        first_level_results.cope_by_state(indind,t,:,f,si)=coapeout;
                                    elseif(strcmp(oat.first_level.cope_type,'cope')),
                                        first_level_results.cope_by_state(indind,t,:,f,si)=copeout;
                                    elseif(strcmp(oat.first_level.cope_type,'acope')),
                                        first_level_results.cope_by_state(indind,t,:,f,si)=abs(copeout);
                                    else
                                        error('Invalid cope_type');
                                    end;

                                    first_level_results.stdcope_by_state(indind,t,:,f,si)=sqrt(varcopeout);
                                    pe(t,:)=peout; % num_timepoints x  num_pes
                                end;

                                %% do bc on the copes
                                if(sum(first_level.bc)>0),
                                    error('No first_level.bc available for continuous data');                            
                                end;
                            end; % freq index                        
                            
                        end;
                        
                    end; % si (state index)

                    %%
                    if(NK>1)
                        %% do fixed effects averaging of states
                        for iContrast = 1:size(first_level_results.cope_by_state,3)
                            % freq loop
                            for iFreq = 1:size(first_level_results.cope_by_state,4)

                                    cope_tmp=permute(first_level_results.cope_by_state(indind,t,iContrast,iFreq,:),[5 1 2 3 4]);
                                    stdcope_tmp=permute(first_level_results.stdcope_by_state(indind,t,iContrast,iFreq,:),[5 1 2 3 4]);

                                    S=stdcope_tmp.^2;

                                    z=ones(size(S,1),1);
                                    [gam, ~, covgam] = flame1(cope_tmp,z,S,1);

                                    first_level_results.cope(indind,t,iContrast,iFreq)=gam;
                                    first_level_results.stdcope(indind,t,iContrast,iFreq)=sqrt(covgam);

                            end % for iFreq = 1:nFreqs
                        end % for iContrast = 1:nContrasts
                    else
                        first_level_results.cope(indind,:,:,:)=first_level_results.cope_by_state(indind,:,:,:);
                        first_level_results.stdcope(indind,:,:,:)=first_level_results.stdcope_by_state(indind,:,:,:);
                    end;
                    
                else
                    first_level_results.cope(indind,:,1,:)=dat(:,:);
                    first_level_results.stdcope(indind,:,1,:)=1;                    
                    first_level_results.cope_by_state(indind,:,1,:)=dat(:,:);
                    first_level_results.stdcope_by_state(indind,:,1,:)=1;
                end;
            
            else,
                if(~isfield(first_level_results,'glm_input_data')),
                    first_level_results.glm_input_times=glm_input_times;
                    first_level_results.glm_input_frequencies=out.tf_freqs;
                    first_level_results.glm_input_data=zeros([Nvoxels_out, size(dat)]);% num_voxels x  num_timepoints x num_contrasts x num_freqs                    
                end;                    
                    
                first_level_results.glm_input_data(indind,:,:,:)=dat;
            end; % if first_level.doGLM                   
        end;
    end; % for indind = 1:Nvoxels_out
    
    if (first_level.doGLM && ~isempty(first_level.connectivity_seed_mni_coord)),
        first_level_results.cope_by_state(first_level_results.seed_index,:,:,:)=1;
        first_level_results.stdcope_by_state(first_level_results.seed_index,:,:,:)=1e16;
        first_level_results.cope(first_level_results.seed_index,:,:,:)=1;
        first_level_results.stdcope(first_level_results.seed_index,:,:,:)=1e16;
    end;
            
    ft_progress('close');
    
    %% save results
    if first_level.doGLM
        first_level_results.x=x;
        first_level_results.contrasts=first_level.contrast;
    end;
    
    first_level_results.source_recon_results_fname=source_recon_results_fname;
    first_level_results.first_level=first_level;
    first_level_results.source_recon=source_recon_results.source_recon;
    first_level_results.level=1;
    %first_level_results.source_recon_results_fname=source_recon_results.fname;
    first_level_results.name=first_level.name;
    first_level_results.recon_method=source_recon_results.recon_method;
    first_level_results.session_name=source_recon_results.session_name;
          
    first_level_results.fname=[ first_level_results.session_name '_' first_level_results.name ];
    disp(['Saving statistics in file ' oat.source_recon.dirname '/' first_level_results.fname]);
    osl_save_oat_results( oat, first_level_results );
    results_fnames{sub}=first_level_results.fname;

    if (strcmp(source_recon_results.recon_method,'none')),
    else
        disp(['To create niftii files from this use a call to osl_save_nii_stats']);
    end;
    
end;

